Finite temperature phase diagram of the classical Kitaev-Heisenberg model 
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We investigate the finite-temperature phase diagram of the classical Kitaev-Heisenberg model on 
the hexagonal lattice. Due to the anisotropy introduced by the Kitaev interaction, the model is 
magnetically ordered at low temperatures for all values of parameters at which the model has a 
discrete symmetry. The ordered phase is stabilized entropically by an order by disorder mechanism 
where thermal fluctuations of classical spins select coUinear magnetic orders in which magnetic 
moments point along one of the cubic directions. We find that there is an intermediate phase 
between the low-temperature ordered phase and the high-temperature disordered phase. We show 
that the intermediate phase is a critical Kosterhtz-Thouless phase exhibiting correlations of the 
order parameter that decay algebraically in space. Using finite size scahng analysis, we determine 
the boundaries of the critical phase with reasonable accuracy. We show that the Kitaev interaction 
plays a crucial role in understanding the finite temperature properties of A2lr03 systems. 
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I. INTRODUCTION 

Recently spin-orbit coupling (SOC) effects have be- 
come a subject of intense research across many different 
disciplines in condensed matter physics. These effects 
are especially pronounced in 4d and 5d transition-metal 
compounds whose significant atomic SOC is due to their 
large atomic weight. 

The strongly entangled, spin-orbital nature of the lo- 
calized states are characterized by an effective angular 
momenta, Jeff. The interactions among the components 
of Jeff are determined by the combination of spin and lat- 
tice symmetries, which creates a possibility for unconven- 
tional anisotropic exchange interactions and, as a result, 
a realization of various novel properties.!"'^ For exam- 
ple, breaking the spin rotation symmetry permits mag- 
netic Hamiltonians that contain terms that are products 
of different components of spin operators. Such terms, 
not allowed in the traditionally studied SU(2) symmetric 
models, introduce a new source of frustration and might 
drive the system towards quantum spin liquid states. 

A prominent model exemplifying these kinds of highly 
anisotropic interactions is the Kitaev model on the hon- 
eycomb lattice J. The ground state of the Kitaev model is 
known exactly; it is a spin liquid characterized by anyonic 
excitations with exotic fractional statistics. Recently, 
Jackeli and Khaliullini suggested that the Kitaev model 
could actually be realized within the honeycomb iridates 
with general formula A2lr03, since the anisotropic part 
of the interactions among Ir-ions has the same form as the 
Kitaev coupling. This suggestion has triggered a lot of 
experimental^"^ and theoretical^l"™ activity in the study 
of Na2lr03 and Li2lr03 compounds. 

In A2lr03 materials, Ir '^ ions are in a low spin 5dP 
configuration with an effective angular momentum of 
Jeff = 1/2 due to strong SOC. The low-energy Hamil- 
tonian describing the interaction between Joff iridium 
moments is called the Kitaev-Heisenberg (KH) model 



which contains both the isotropic antiferromagnetic (AF) 
Heisenberg and the anisotropic Kitaev interactions. 

The Kitaev exchange interaction is generated through 
the 90° Ir-O-Ir hopping path between Jeff Kramers dou- 
blets and is non-zero only in the presence of finite Hund's 
coupling. At the same time, the isotropic AF Heisenberg 
interaction is suppressed; the only non- vanishing contri- 
bution is from the superexchange interaction arising from 
the direct overlap of the Ir 5d orbitals. The isotropic ex- 
change via the 90° Ir-O-Ir is canceled due to a destructive 
interference among multiple superexchange paths. 

The ground state of the KH model can be a spin liq- 
uid despite the presence of the direct isotropic exchange. 
This is because the Kitaev spin liquid is rather stable 
with respect to the Heisenberg interaction.2ii^; it remains 
the ground state of the KH model for a wide range of the 
Heisenberg interaction. Even so, the Kitaev spin liquid 
in honeycomb iridates has not been observed yet; in both 
Na2lr03 and Li2lr03, magnetic order has been observed 
at low temperatures!^"— Moreover for Na2lr03, the KH 
model in its original form does not appear to be sufficient 
to account for neither the zigzag magnetic order nor the 
spectrum of magnetic excitations measured in neutron 
scattering experiments.—"— 

To explain the experimental observations in Na2lr03, 
three different modifications of the super-exchange model 
describing the interactions between localized magnetic 
moments have been proposedJi"^^ In the first ap- 
proach, ''^^ it was shown that the zigzag magnetic order 
may be stabilized within the KH model by including sub- 
stantial second and third neighbour antiferromagnetic 
interactions. The second approacbi^, extends the KH 
model to its full parameter space by including additional 
hopping processes based on the t2g — eg hopping along 
the 90° Ir-O-Ir paths. The main difference between these 
two approaches is that the role of the Kitaev interac- 
tion is minor in the first approach while it still plays the 
dominant role in the second one. The third approach 
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FIG. 1: Four possible magnetic configurations: (a) FM or- 
dering; (b) two-sublattice, AF Neel order; (c) stripy order; 
(d) zigzag order. Open and filled circles correspond to up 
and down spins. 



assumes that Na2lr03 is significantly distorted from tlie 
ideal structure; liere, tlie trigonal distortion is the dom- 
inant interaction while the SOC is subdominant. This 
third model actually does not contain the Kitaev term 
at all, but it does retain the zigzag ground state as one 
of its magnetic ground states.™ Thus, the importance of 
the Kitaev term in the low-energy physics of the A2lr03 
compounds remains an open question. 

In this paper, we continue to examine the role of the 
Kitaev interaction in the magnetic properties of A2lr03 
systems using the KH model defined in both the orig- 
inal^ and the extendedi^ parameter space. We argue 
that the Kitaev interaction's main effect is the reduc- 
tion of the symmetry of the system from the continu- 
ous SU(2) symmetry to a discrete Zg symmetry. The 
latter is crucial for the finite temperature properties of 
quasi two-dimensional (2D) iridates since by the Mcrmin- 
Wagner theorem^S the 2D magnetic systems with contin- 
uous symmetry do not exhibit long-range magnetic order 
at any finite temperature. However, there is no such con- 
straint on spin systems with discrete symmetry. Thus, we 
argue that the Kitaev term is responsible for the presence 
of the long range magnetic order at temperatures which 
are significantly higher than those that one can expect 
from a weak inter-layer coupling making A2lr03 com- 
pounds weakly three-dimensional. 

We also show that the finite temperature properties 
of the KH model are similar to those of the six-state 
clock model. ^^'^^ The KH model undergoes two continu- 
ous phase transitions as a function of temperature. This 
gives rise to three different phases: a low-T ordered phase 
with a spontaneously broken Zg symmetry, an intermedi- 
ate critical phase, and a high-T disordered phase. Finite 
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Classical energy as a function of a: in tiie FM phase 
3 — 5a, dot-dashed green line), in the zigzag phase 
-3a + 1, dotted black line), in the Neel phase (S^ = 



5a — 3 solid blue line), and in the stripy phase (E^i — 
dashed red line). 
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size scaling analysis of our Monte Carlo (MC) simula- 
tions confirmed that the intermediate phase is a critical 
Kosterlitz-Thouless (KT) phase with floating exponents 
and algebraic correlations. 

The rest of the paper is organized as follows. We start 
our discussion with Sec. HI] which contains a brief sum- 
mary of known facts about the classical KH model that 
will be used in the rest of the paper. In Sec. IIIIl we dis- 
cuss in detail the results of our numerical simulations. We 
analyze the discreteness of the low-temperature phase, 
and then we discuss the criticality of the intermediate 
phase and show how the finite-size scaling analysis allows 
us to determine the boundary of the critical phase. We 
show that the phase transitions are driven by the order by 
disorder mechanism, in which the thermal fluctuations of 
classical spins select collinear spin configurations where 
the magnetic moments all point along one of the cubic 
directions. We also demonstrate that at particular points 
of the phase diagram, for which the continuous symmetry 
is preserved, the ordered state is destroyed at any non- 
zero temperature. The section ends with a discussion of 
the finite temperature phase diagram of the classical KH 
model. 

In Sec. lIVl we analyze the finite temperature properties 
of the KH model in its extended parameter space. We 
show that for the parameters relevant to the Na2lr03 and 
Li2lr03 compounds, the model exhibits a zigzag mag- 
netic order in agreement with experimental findings. In 
Sec|Vl we show that there is a significant difference be- 
tween the KH model and the Heisenberg model with a 
cubic anisotropy. While the later model also exhibits two 
phase transitions, its intermediate phase is not critical; 
it is nematic like. Using finite-size scaling analysis, we 
show that in this case the phase transitions are of the 
three-states Potts and Ising universality classes. Sec. IVII 
contains a summary of the obtained results and conclu- 



(111) 




(b) 



(111) 



(c) 



(111) 



(d) 



(111) 



^.hi 







FIG. 3: A low-T distribution of the projections of the vector order parameter on the (111) plane, a) Neel order parameter for 
a = 0; b) Neel order parameter for a = 0.25; c) stripy order parameter for a — 0.5; d) stripy order parameter for a = 0.75. 



II. THE CLASSICAL KH MODEL 



The classical version of the KH model is 
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Following the notation of Ref. 13, we denote the effective 
degrees of freedom of Ir*+ ions as S, and take the quan- 
tization axes along the cubic axes of the IrOe octahedra. 
J = x,y, z denotes the three bonds of the honeycomb lat- 
tice. Two couplings, Jh and Jk, are opposite in sign; the 
AF isotropic exchange, Jh, and the ferromagnetic (FM) 
anisotropic exchange, Jk. 

The exchange constants corresponding to the Kitaev 
and the Heisenberg interactions in the KH model ([T]) can 
be conveniently described by one parameter, a, such that 
Jk = '2a and Jh = I — a. The model then reads as 
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Let us briefly describe possible magnetically ordered 
states of model ^. For this we need to introduce four 
sublattices: A, B, C, and D. For < a < 1/3, the clas- 
sical ground state is the simple two-sublattice AF Neel 
order (Fig.l (b)) characterized by the order parameter 
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Af denotes the number of sites. The stripy AF order 
(Fig.l (c)) describes the classical ground state of the 
model for 1/3 < a < 1 and is given by 
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For a ~ 1.0, the stripy AF state is classically degener- 
ate with other magnetically ordered states (see Figl2]). 
For example, it is degenerate with FM order (Fig.l (a)) 
described by the total magnetization 
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and with zigzag AF spin order (Fig.l (d)) described by 
Z = -r-: ^(SiA + Sis - Sic - Sio) ■ (6) 



The classical degeneracy of the point corresponding to 
the classical Kitaev model is macroscopic and is known 
exactly; asymptotically, it has (1.662)^ spin configura- 
tions. This was computed by Baskaran et aP^ by map- 
ping ordered states of the classical Kitaev model to a 
certain dimer covering of the honeycomb lattice whose 
total number is knowni^ 

Because of the presence of the anisotropic Kitaev in- 
teraction term, the KH model has discrete spin-rotation 
symmetry for all of a except at two points, a — and 
a — 0.5. Ata = 0, the KH model reduces to the anti- 
ferromagnetic Heisenberg model with continuous SU(2) 
symmetry. At a = 0.5, the stripy phase becomes an ex- 
act ground state of the model; it corresponds to a fully 
polarized FM state in a rotated basis. This basis can be 
seen by fixing the spin's direction on sublattice A and ro- 
tating the spins on sublattices B, C, and D by the angle 
TT about the x, y, and z axis, respectively (see Fig.[TJc))j^ 
It is evident that the FM state has true SU(2) symme- 
try. Away from these special points, the symmetry of the 
KH model is discrete; it combines the cubic symmetry of 
both the spin and the lattice space. 



III. FINITE TEMPERATURE PHASE 

DIAGRAM AND CRITICAL PROPERTIES OF 

THE CLASSICAL KH MODEL 

In this section we study the behavior of the classical 
KH model ([2]) at finite temperature using MC simula- 
tions based on the standard Metropolis algorithm. In 
our simulations, we treat the spins as three-dimensional 
(3D) vectors, S^ = {Sf,Sf,S^), of unit magnitude with 
{Sfy + {Sf)^ + {S^f = 1. At each temperature, more 
than 10^ MC sweeps were performed. Of these, 10^ were 
used to equilibrate the system, and afterwards only 1 out 
of every 5 sweeps was used to calculate the averages of 
physical quantities. We present all energies in units of 
Jh and assume fc^ = 1- The simulations were performed 
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FIG. 4: Snapshots of spin configurations in the low-T phase computed for a) q = 0, b) q = 0.25 c) a = 0.5 and d) a — 0.75. 
Each color corresponds to one of the two or four different sublattices describing the Neel or the stripy order parameters, 
respectively. 



on different systems with a total number of sites equal 
to 2 * L * L. The systems are spanned by the primi- 
tive vectors of a triangular lattice ai = (1/2, V^/2) and 
^2 = (1,0) with a 2-point basis using periodic boundary 
conditions. 

In our simulations we compute the following observ- 
ables: four different order parameters O = {N, S, M, Z} 
defined in Eqs.(l3][6]), corresponding susceptibilities 



Xo=AA((02)-(0)2)/T, 
the Binder's cumulants 

and the specific heat 

C = {{E^) - {Ef)/NT^ . 

A. Low-temperature ordered phase 
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At low temperatures, the KH model magnetically or- 
ders in either a Neel state, N, or in a stripy state, S, 
depending on the relative strengths of the Kitaev and 
the Heisenberg interactions. The presence of long range 
order at finite temperatures requires a discreteness of the 
order parameter, and means that the direction of the or- 
der parameter must also be selected. This, however, does 
not happen on the level of non-interacting spin waves. 
For both the Neel state and the stripy state, the linear 
spin-wave spectrum has a quasi-Goldstone mode at the 
ordering vector simulating the spontaneous breaking of 
the continuous symmetry.^ The discrete symmetry of the 
order parameter appears only due to the contribution 
of higher order anharmonic modes of spin fluctuations 
which remove the accidental degeneracy of the classical 
ground states and lower the energy of the states whose or- 
der parameter points along cubic directions.^ This gives 



a six-fold degeneracy of the order parameter manifold 
which corresponds to the six degenerate minima in the 
free energy. This is a well known order by disorder mech- 
anism in which spin fluctuations (quantum or thermal) 
remove the accidental degeneracy and select the true or- 
dered state. 

The discreteness of the order parameter at all but spe- 
cial points can be revealed from the histogram method 
in which we record two-dimensional distributions of the 
projections of the vector order parameter on the (111) 
plane. We use the projection of the order parameter on 
the (111) plane because it preserves the cubic symmetry 
of the model and shows that its degeneracy is related to 
the orientation of the order parameter with respect to 
the directions in a cubic crystal. In Fig. [3] (a), (b), (c), 
and (d) we present, respectively, the low-T distributions 
of the projection of the vector order parameter on the 
(111) plane for four different values of a: 0, 0.25, 0.5, 
and 0.75. The six-fold peak structure is observed in the 
distribution function of the Neel and the stripy order pa- 
rameters for a = 0.25 (Fig. [3] (b)) and a = 0.75 (Fig. [3] 
(d)). The peaks correspond to the ordered spin config- 
urations with an order parameter pointing along one of 
the cubic axes. In Figs. 0] (b) and (d) we present snap- 
shots of the spin configurations which contribute to the 
histograms in Figs. |3] (b) and (d). In the snapshot of 
the spin configuration shown in Fig.|31Jb), the spins that 
are pointing along the y— direction are antiparallel on 
the two sublattices. This state corresponds to the Neel 
phase with the order parameter directed along y— axis. 
In the spin configuration of Fig. 2] (d), spins are point- 
ing parallel to the z— direction on two sublattices and 
antiparallel on another two. This state corresponds to 
the stripy phase with the order parameter directed along 
z— axis. 

The projections and snapshots of the order parameter 
for the two special points a = and a — 0.5 are pre- 
sented in Figs. [3] (a) and (c) and Figs. S] (a) and (c). 
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FIG. 5: Histograms of the order parameter 7njv(s)j obtained for the system with 2*84*84 spins in the ordered phase, (a) and 
(e), in the intermediate phase, (b)-(c) and (f)-(g), and in the disordered phase, (d) and (h). Histograms (a)-(d) are computed 
for a = 0.25, and (e)-(h) are for a — 0.75. The histograms are presented on the complex plane (Re |mjv(s)|, Im |7n]v(s)l)- 



respectively. These are the points with continuous sym- 
metry. As expected, we see that the projections of the 
order parameter on the (111) plane (|3] (a) and (c)) are 
more or less equally distributed along a circle. Note that 
the multiple " quasi" -peaks appear because the computa- 
tions are performed on a finite size system. The snapshot 
for Q = (HI (a)) shows that a certain direction is chosen 
in this particular spin configuration but it is not along 
one of the cubic axes. Because this direction varies from 
configuration to configuration and also because the spins 
belonging to the same sublattice can be orientated paral- 
lel or antiparallel to this direction, there is no long range 
order for a = 0.0. A similar situation is observed at 
a = 0.5. Here, spins at each sublattice point along a cer- 
tain direction which is not chosen to be along any cubic 
axes. There is no long range magnetic at a = 0.5 as well. 



Finally, there is no additional degeneracy of the or- 
der parameter related to its real-space structure. Both 
the stripy and the zigzag order can have three different 
configurations of bond orderings due to the 120° rota- 
tional symmetry of the honeycomb lattice. The absence 
of the additional degeneracy can be understood from the 
fact that, due to SOC, the symmetry transformations 
act simultaneously on spins and on the lattice. Once the 
pattern of the bond ordering is chosen, the direction of 
the order parameter in spin space is also chosen. For 
example, in the vertical stripy state shown in Fig. [Ijc) , 
the spins are directed along the z axis just like in the 
snapshot shown in Fig.|3](c). 



B. The critical nature of the intermediate phase. 
Finite size scaling analysis. 



Regardless of the specific kind of magnetic order, the 
low-T ordered phase is separated from the high-T para- 
magnetic phase by the intermediate phase. In our re- 
cent study, we have shown that the intermediate phase 
is a critical phase with two finite-temperature bound- 
aries that correspond to Berezinskii-Kosterlitz-Thouless 
(BKT) phase transitions.''^^ 

To describe the finite temperature properties of the 
KH model, we use a projection of the vector order pa- 
rameter on the (111) plane. The vector is character- 
ized by both the absolute value and the azimuthal angle 
of |?7iAr(s)|. This projection is equivalent to the Zg or- 
der parameter, and we write it in a complex form via 
™Af(s) = J2i=i I'^i.AfCSjIe**'. We chose the phase 9 such 
that the minimal-energy states of the order parameter, 
which point along the cubic axes, will be labeled by the 
values 9i = ttoj/S, ni = 0,..5. Observing the critical 
phase proved to be challenging for several reasons. First, 
within the vicinity of the BKT transition, the critical be- 
havior gives rise to a very slow, logarithmic convergence 
to the thermodynamic limit. We thus had to perform the 
finite size scaling analysis of our simulations on rather 
large systems with L = 84,96,108,120,144,168,204. 
Second, contrary to the Ising-like spin systems in which 
all previous studies of the six-state clock model's critical 
phase have been performed, the magnetic degrees of free- 
dom in the classical KH model are 3D Heisenberg spins 
which are affected by thermal fluctuations more strongly. 
As a result, a large number of sweeps is needed to aver- 
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FIG. 6: The log-log plots of the order parameter mjsii^s) 
as a function of system size L at various temperatures. The 
sohd curves indicate the linear behavior that corresponds to a 
power-law dependence, mjv{s) ~ I/"''", corresponding to the 
intermediate critical phase. The dashed curves show deviation 
away from the linear behavior outside the critical phase. 
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FIG. 7: The Binder cumulant as a function of temperature 
for (a) a = 0.25 and (b) a — 0.75. From the crossing points of 



different Binder's curves, we estimate Tc, 



: 0.152 and Tc, 



0.124 for a = 0.25 and a — 0.75, respectively. 



age out the thermal fluctuations and capture the critical 
behavior. 

A clear indication of the three-phase structure emerges 
from the inspection of the histogram plot of the com- 
plex magnetization mpf(^s^ at different temperatures (see, 
Fig|S]). As we move from low to high temperature, we 
observe a transition from an ordered phase (6 isolated 
spots) through an intermediate critical phase (ring dis- 
tribution) up to the disordered phase (uniform distribu- 
tion around zero). The critical phase has an emergent, 
continuous U{1) symmetry)^ which is reminiscent of the 
intermediate phase of the six-state clock model.?J- At the 



boundaries, Tc^ and Tc2, and inside the critical phase, 
the order parameter exhibits a power law dependence on 
system size of the form, tojv(s) ^ L^^^^. From renormal- 
ization group analysis;^ it is known that the lower Tci 
and the upper Tc2 transitions in the six state clock model 
are characterized by the critical exponents 771 = 1/9 and 
V2 = 1/4, respectively. 

Fig. [S] shows the log-log plots of the order parameter 
TOjv and ms as a function of system size L for different 
temperatures. In accordance with the power law behav- 
ior of the order parameter, the data points of the log-log 
plot show a linear behavior inside the critical phase. For 
a = 0.25 (Fig.|n](a)), the log-log plots of the order param- 
eter as function of L show a linear behavior in the tem- 
perature interval between Tc-^ ~ 0.152 and Tc^ — 0.162, 
with the boundaries characterized by the critical expo- 
nents 771 = 0.13 and 772 — 0.21. For a — 0.75, we detected 
the critical phase in the temperature interval between 
Tci ~ 0.125 (?7i = 0.13) and T^^ ^ 0.127 (7/2 = 0.22). For 
both a — 0.25 and a = 0.75, our estimates for the criti- 
cal exponents of the boundaries are in a good agreement 
with theoretical values for the six state clock model. 

We can verify our estimates for the BKT transition 
temperatures using the Binder cumulants method and 
further finite-size scaling analysis.?- The Binder cumu- 
lants Bm^ns) 1 defined in Eq.®, have a scaling dimension 
of zero; thus the crossing point of the cumulants for dif- 
ferent lattice sizes provides a reliable estimate for the 
value of the critical temperature T,.^ at which the long 
range order is destroyed. The results for the Binder cu- 



mulants B 



™JV(S) 



are presented in Fig. [T] The crossing 
points for a = 0.25 and a = 0.75 are at Tc;^ = 0.153 and 
at Tci = 0.125, respectively, and are in a good agreement 
with estimates obtained from the log-log plots in Fig. |6l 
Another way to check the nature of the critical phase 
boundaries is by using the data collapse method based on 
the finite size scaling arguments. In the BKT transition, 
the order parameter and its susceptibility exhibit power 
law behavior, mf^(^g'j ^ ^-v/^ ^nd Xn{S) ^ S,^^^ ^ while 
the correlation length near the critical temperature, Tc, 
diverges as ^ ~ exp(ai^"), where a is a non-universal 
constant and t — \T — Tc\/Tc is the reduced tempera- 
ture."^^ The finite size scaling analysis is based on the 
assumption that the singular part of the free energy is a 
homogeneous function of system size L and correlation 
length ^, and depends only on their ratio L/^. Based 
on this assumption, the finite size scaling behavior of the 
order parameter and the susceptibility have the following 
functional forms 
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(10) 



Amjv(S) ~ ^ ""iJV(S) 



where scaling constants h — 77/2 and c = 2 — 77, and 
M^ig\ and S^^^ are unknown universal functions. We 
plot the variation in temperature of the order parameter 
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FIG. 8: Thermodynamic quantities and finite size scaling analysis, (a) and (b): The order parameters tun and ms as functions 
of T for a = 0.25 and a — 0.75, respectively, (c) and (d): The susceptibility, X™iv(s) (defined in Eq.lJT])) plotted as a function of 
T for a — 0.25 and a = 0.75, respectively, (e) and (f): Finite-size scaling of the order parameter data in the low-temperature 
region (T < T^) for a — 0.25 and a = 0.75, respectively, (h) and (g): Finite-size scaling of the high-temperature susceptibility 
data ( T > Tc2) for a — 0.25 and a = 0.75, respectively. 



in FiglH] (a) and (b) and the susceptibility vs temperature 
in FiglH] (c) and (d). In FigE] (e) and (f) and FigE] (h) 
and (g), we show the scaling plots for m]^ig\L^ as a func- 
tion of L^i cxp(a/^(Tci -T)/TcJ and Xm„(s)i"'' as a 
function of L~-^ exp(a/Y/(T^^T^JJ/7^), respectively. 

The finite-size effects in the behavior of the order pa- 
rameter and its susceptibility are striking: we see both a 
significant finite-size tail of mpf/g\ into intermediate re- 
gion for T > Tc-^ and a strong dependence of the position 
and height of the susceptibility peak on the system size 
L. This is a consequence of the infinite correlation length 
in the intermediate critical region, which in the finite size 
systems look like the quasi long range order. Neverthe- 
less, the data points for different system sizes plotted 
in the scaled form collapse reasonably well on universal 
curves that correspond to the universal functions Af7v(s) 
and 2mjv(s) ■ ^^^ best data collapse was obtained for the 
following scaling parameters: a — 1.9, b = 0.056 and 
c = 1.45 and transition temperatures Tc^ = 0.153 and 
Tc, = 0.1615 for a = 0.25; a = 1.55, b = 0.056 and 
c — 1.55 and transition temperatures Tc-^ = 0.125 and 
Tc^ = 0.127 for a = 0.75. The values obtained for b and 
c give the following critical exponents for the lower and 
the upper boundary of the critical phase: ri{Tc^) — 0.11, 
rj{Tc^) = 0.275, and r]{Tc,) = 0.11, 77(TcJ = 0.225 for 
a = 0.25 and a = 0.75, respectively. 

To summarize the discussion of the numerical data pre- 
sented above, we can say that we definitely observe the 
critical intermediate phase in the classical KH model, 
which, however, slightly deviates from the standard BKT 
criticality. We believe that the imperfect data collapse to 



be caused by the presence of three dimensional spin fluc- 
tuations of Heisenberg spins and finite sized system sizes. 
Nevertheless, the smallness of the discrepancy and over- 
all high quality of the scaling indicate that this effect is 
subdominant in the temperature range where the critical 
behaviour occurs. 



C. The finite temperature phase diagram 

In Figl9] we present the finite temperature phase dia- 
gram of the KH model ([2]) obtained by tracking the de- 
pendance of the transition temperatures Tc-^ and Tc^ on 
the strength of the KH model parameter a. As discussed 
in detail in the previous section, in order to determine 
critical temperatures we have used the Binder's cumu- 
lant method and finite size scaling arguments. 

As expected, only two low-T magnetically ordered 
phases are present on the phase diagram - the Neel phase 
at a < 1/3 and the stripy AFM phase at 1/3 < a < 1. 
The vertical line that separates the Neel and the stripy 
AFM phases corresponds to a first-order phase transi- 
tion. Our numerical results show that the critical phase 
is present for all values of a except the special points 
a — Q, 0.5, and 1 at which the model has continuous 
symmetry, and the ordered state is destroyed at any non- 
zero temperature in accordance with Mermin- Wagner 
theorem.Sfl: Note that the way the width of the critical 
phase disappears in the vicinity of a = and a = 0.5 
is different than for a = 1. In Fig 1^1 we can see that 
the critical phase narrows very rapidly in the vicinity of 
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FIG. 9: Phase diagram of the classical KH model ([2]). The 
regions designated by, "N", and "S", are the regions of the 
phase diagram where the Neel, and stripy order persist. The 
vertical line between the Neel and the stripy phase at a = 
1/3 denotes the first order phase transition. Each red circle 
designates a lower critical temperature Tc;^ for the value of 
a for which the model was explicitly simulated; the red line 
extrapolates between these points. For each a, the critical 
value of Tc-^ was determined through the crossings of Binder's 
cumulant curves. The critical phase is shown in blue. The 
upper boundary of the blue critical region was determined by 
finding the temperature Tc2 , for which the value of the critical 
exponent rj exceeds 0.25. The blue line extrapolates between 
each computed point. 



a = 0,0.5. This behavior is related to the fact that at 
these particular points the anisotropy either does not ex- 
ist as at a = or cancels out as at a = 0.5. In the 
vicinity of these points, the local minima of the free en- 
ergy that corresponds to different directions of the order 
parameter are still separated by finite energy barriers. 
These barriers allow finite temperature ordering through 
the order by disorder mechanism. At the same time, the 
width of the critical phase slowly decreases when we ap- 
proach the Kitaev limit and is extremely small in the 
vicinity of a = 1. In this limit, the difference between 
the classical energies of the stripy, zigzag and ferromag- 
netic phases is decreasing with the critical temperature 
tending to zero. At a = 1, all of the phases are degener- 
ate (see also Fig. [5] and the discussion above), thus, there 
is no order- by-disorder in this limit. 



IV. FINITE TEMPERATURE PHASE 

DIAGRAM OF THE EXTENDED CLASSICAL 

KH MODEL 

Here we explore the finite temperature phase dia- 
gram of the KH model extended to its full parameter 
space. This extension, consisting in taking into account 
all super-exchange processes leading to the coupling be- 
tween Ir ions, gives a possibility to capture the zigzag 
magnetic order observed in the A2lr03 compoundsjii^ 

There are three physically distinct processes that de- 
termine the ratio between the Kitaev and the isotropic 



Heisenberg terms in the original KH model. Two of the 
processes involve the virtual hopping of the t2g electrons 
through the two nearest oxygen ions. As it was shown 
by Jackeli and Khaliulin- , the processes via the upper 
and lower oxygen interfere destructively and the isotropic 
part of the Hamiltonian exactly vanishes. Exchange cou- 
plings of neighboring Kramers states on Iridium ions ap- 
pear due to either the multiplet structure of the excited 
levels on a Ir ion induced by Hund's coupling, or due to 
Coulomb interactions. The third process involves a di- 
rect hopping between NN t2g orbitals which gives a finite 
Heisenberg term in the KH model. All of these processes 
involve only t2g electrons. However, there is another pos- 
sible process)^ the intersite t2g ^ eg hopping along the 
90° Ir-O-Ir paths. This is the dominant pathway in a 90° 
geometry since it involves strong tpda overlap between the 
p— orbitals on oxygen and the Cg orbitals on iridium ions. 
Remarkably, these hopping processes also introduce the 
Kitaev interaction, but with a different sign. 

Following Refll9l, we consider the most general exten- 
sion of the KH model, in which both Kitaev and Heisen- 
berg interactions can change sign and be both FM and 
AFM. The Hamiltonian of the extended model is 



n = A(2sin( 



' L SIS] 






fll) 



where the relative strengths of the effective interaction 
between Ir magnetic moments are described by the phase 
angle, defined in the interval (0,27r), and the overall 
energy scale, A — \/ J\ -|- J^. The phase space of the 
original KH model is covered by the values of </) from 
37r/2 to 27r. 

Using MC simulations, we investigate the finite tem- 
perature properties of this extended model (llip . and 
present the finite temperature phase diagram in Fig llOl 
In its full parameter space, the extended KH model ac- 
commodates four different classical phases. In addition 
to the Neel and the stripy phases present in the phase 
diagram (FiglH]), the FM and the zigzag phases appear 
in Fig llOl in a wide range of values of (/>. The previously 
missing zigzag magnetic order is found to occupy almost 
a quarter of the phase space of the extended model. Our 
findings are in agreement with the ground state phase 
diagram of the quantum model ([TT|) obtained by ex- 
act diagonalization..^^ However, as we already discussed 
above, the classical model does not support the spin liq- 
uid phases present in the quantum model. 

The finite temperature properties of the extended 
model are very similar to the finite temperature prop- 
erties of the original KH model. Both the transition be- 
tween the zigzag and the FM order, and the transition 
between the Neel and the stripy are first order phase 
transitions. The intermediate phase is present for the 
whole parameter space of the extended model except for 
special points at which the continuous symmetry is re- 
stored. The similarities between finite temperature prop- 
erties of the original and the extended KH model are 
not surprising. The above mentioned 4-sublattice spin 



is defined by the Hamiltonian 
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FIG. 10: Phase diagram of the extended classical KH model 
(fTTjI . The regions designated by, "N", "S", "Z", and "F" 
are the regions of the phase diagram where the Neel, stripy, 
zigzag, and FM order persist. Each red circle designates a 
lower critical temperature T^ for the value of for which 
the model was expUcitly simulated. For each value of (f}, the 
critical value of Tc^ was determined through the crossings of 
Binder's cumulant curves; the red line extrapolates between 
these points. The critical phase is shown in blue. The upper 
boundary of the blue critical region was determined by find- 
ing the temperature Tc2 , for which the value of the critical 
exponent rj exceeds 0.25. The blue line extrapolates between 
each computed point. 



transformation^^ permits a mapping of the AF stripy 
phase to the FM phase and of the AF Neel to the zigzag 
phase. This implies that these phases have identical dy- 
namics at both zero temperature^^ and at finite tempera- 
tures, which is in agreement with our numerical findings. 
According to the fit of the uniform magnetic suscepti- 
bility presented in ReiM, ~ 111° ± 2° ~ 0.627r± O.OItt 
for NailrOa and (f> ~ 124° ± 4° ~ 0.697r ± O.OStt for 
Li2lr03. Using these values, we get the following esti- 
mates for the Neel temperatures: Tc-^ ~ 0.16 which is 
equivalent to about 17.7 K for Na2lr03, and T^^ ~ 0.19 
which is equivalent to about 21 K for Li2lr03. Both 
values are close to the experimental value Tjv ~ 15 K ob- 
tained for both Na2lr03 and Li2lr03 compounds. ^'^ Our 
estimates for the upper boundary of the critical phase 
is Tc, = 0.18 for Na2lr03 and T^, = 0.22 for Li2lr03 
which is equivalent to about 20 K and 24.5 K, respec- 
tively. This gives the estimates for the width of the criti- 
cal phase about 2.3 K and 3.5 K for Na2lr03 and Li2lr03, 
respectively. 



CLASSICAL HONEYCOMB HEISENBERG 
ANTIFERROMAGNET WITH A CUBIC 
ANISOTROPY 



To complete this study, we also investigate the nature 
of the finite-temperature phase transitions of the hon- 
eycomb Heisenberg antiferromagnet in which the cubic 
anisotropy is included explicitly. In this case, the model 






E 



{{s:r+{sfr+(sn') (12) 



where D denotes the strength of the cubic anisotropy. 
The positive sign of the anisotropy, D > 0, determines 
that spins tend to align along the cubic axes, and not in 
the diagonal directions of the lattice which would corre- 
spond to I? < 0. Because of the cubic anisotropy, the 
model (|12|) shows no full rotational symmetry. The spin 
projections in the (111) plane again will favor six direc- 
tions as in the case of the KH model. Here we explore 
by MC simulations whether or not the finite temperature 
properties of models ([T]) and ([T2|) are similar. 

Let us discuss the thermodynamic properties of the 
model (mi). At T = and in the absence of the 
anisotropy, the ground state of the model (fT2)) is deter- 
mined by the Neel order parameter N. In the presence 
of the anisotropy D, the ordered state will survive un- 
til some critical temperature. At high temperature, the 
thermal fluctuations of spins will destroy the magnetic 
order and the system will be in the disordered paramag- 
netic state. As in the KH model, in order to reach the 
ordered state from high temperature, the system has to 
break the discrete Zg symmetry of the Heisenberg honey- 
comb model with cubic anisotropy. In the case of the KH 
model, the Zg symmetry was broken through the critical 
intermediate phase. In our simulations of (jl2[) . we did 
not observe any critical phase; instead, the Zg symmetry 
is broken in the following two steps. After lowering the 
temperature, the Z3 symmetry is broken first at T = Tc^'- 
the system still remains paramagnetic since the expecta- 
tion value of (N) = 0. In the second step, the remaining 
Z2 symmetry is broken at the temperature T = T^ < Tc^ 
and the system acquires long-range magnetic order. 

The numerical results computed for the value of the 
cubic anisotropy D = 0.1 are presented in Figs. [11] (a)- 
(d). They display the data for the temperature depen- 
dence of the staggered magnetization N and the cubic 
parameter bjq and corresponding susceptibilities. In our 
simulations, the cubic order parameter 6jv is expressed 
as the expectation value of a doublet given by 



bNi = {Nl + Nl - 2Nl)/./& , 
bN2 - [Nl - N^)/V2 , 



(13) 
(14) 



where N^ , Ny , N^ are the components of the Neel order 
parameter N. At high temperature, the cubic symmetry 
is not broken and N^ = Ny = N^ and (bNi) = {bN^) = 0- 
As we can see from the Figs. [11] both the order parame- 
ters and their susceptibilities indicate a continuous phase 
transition at different temperatures Tc^ and Tc^ for N 
and &Ar, respectively. At the higher temperature, T^^, 
the cubic symmetry is spontaneously broken and one of 
the cubic axes is selected, say z. Then {bj^i) = — ■\/2/3, 
(67V2) — and the cubic parameter (brvf) acquires a fi- 
nite value. Nevertheless, there is still no long range order 
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FIG. 11: The magnetization A'^ (a) and the cubic order parameter bn (b) as functions of T. The susceptibihties XJV (c) and 
Xbjv ('i) ^^ functions of T. (e) and (f): Finite-size scaling of the order parameters A^ and 6jv in the low-temperature region 
T < Tci- (h) and (g): Finite-size scaling of the high-temperature susceptibilities xn and Xbjv in the high-temperature region 

0.1. 



as the time reversal symmetry remains unbroken. Thus, 
we can say that the intermediate phase is nematic-like. 
At the lower temperature, Tci , the time reversal sym- 
metry is also broken and the system acquires long-range 
magnetic order characterized be non-zero (N). We can 
estimate Tci and T^^ using the associated Binder's cumu- 
lants, Bn and B^^ , whose temperature dependencies are 
presented in Figs[T2] (a)-(b). From their crossing points 
we estimated the transition temperatures to be equal to 
Tci ^ 0.314 and T^^ = 0.32. 

Next, in order to obtain the critical exponents charac- 
terizing these two phase transitions we perform a finite- 
size scaling analysis. As the Zg symmetry of the model 
(|12p is reduced to Z2 symmetry in the intermediate 
phase, the high-T transition is expected to be in the 
universality class of 3-states Potts model and the low- 
T transition to belong to the Ising universality class. 

First, using the scaling relation for the Binder cumu- 
lants, Bn and Bh^^, we obtained the correlation length 
exponents i^i and 1^2 describing the divergence of the cor- 
relation length, ^ ^ \T ~ Tc\''^'^''\ close to the critical 
points Tc — Tc^ and Tc — Tc^, respectively. In accor- 
dance with the theory prediction, the best data collapse 
is obtained for z^i = 0.83 and 1^2 = 1-0 which are the 
critical exponents of the three-states Potts and 2D Ising 
model, respectively. 

Second, having found the critical exponents I'l and 1^2, 
we can estimate the critical exponents /3 and 7 by per- 
forming a scaling fit of the order parameters and suscep- 
tibilities (Figs. [TT] (e)-(g)). Near the critical tempera- 
tures, Tci and Tc^, the cubic order parameter and stag- 
gered magnetization are expected to satisfy the scaling 



relation 5jv ^ L^^^/"^ and TV ~ L^^'il^'^^ respectively. 
The scaling laws for their susceptibilities are given by 
X6„ ~ LTi/''! and xn ^ L^^/"-'. The best scaUng is ob- 
tained for the exponents /3i = 1/9 and 71 = 13/9 - the 
exact scaling coefficients for the Pott's transition, and 
P2 — 1/8 and 72 — 7/4 - the exact scaling coefficients for 
the Ising transition. 

To conclude this section, in Fig. 13 we present a finite 
temperature phase diagram for the model p^ . We see 
that the width of the intermediate phase decreases with 
the increase of the value of the cubic anisotropy parame- 
ter D. At around D ~ 0.5, two transitions become indis- 
tinguishable and the intermediate phase collapses. The 
transition between the low-T ordered magnetic phase and 
high-T disordered phase becomes first order. 



VI. CONCLUSION 

In this paper we studied finite temperature properties 
of the classical two-dimensional KH model and computed 
the phase diagram of this model in its full parameter 
space. We started by analyzing the lowest-energy mag- 
netic configuration and then found that all of the mag- 
netic phases are accompanied by an accidental continu- 
ous rotational degeneracy which does not correspond to 
any symmetry of the Hamiltonian. This pseudo degener- 
acy is lifted by thermal fluctuations of spins giving rise 
to an ordering at low temperature as observed by our 
MC simulations. Specifically we determined that the low 
temperature phase is magnetically ordered at all values 
of parameters at which the model has discrete symme- 
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FIG. 12: (a) and (b): The Binder cumulants, Bn = 1 - 
{N^)/3{N'^f and Bt^ = 1 - {b%)/3{b%f, as functions of 
T. (c) and (d): Finite-size scaling of the Binder cumulants 
Bn and Btj^ . The anisotropy constant is considered to be 
D = 0.1. 
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FIG. 13: Finite temperature phase diagram for the classi- 
cal honeycomb Heisenberg antiferromagnet with the cubic 
anisotropy <\l'2\i . The blue line denotes the high-T three state 
Pott's transition, and the red line denotes the low-T Ising 
transition. 



try. The direction of the order parameter is chosen by 
thermal fluctuations of the spins through the order by 
disorder mechanism. 

From numerical MC simulations, we have verified that 
the classical KH model effectively behaves like the six- 
state clock model. At finite temperatures, the model 
exhibits two phase transitions and an intermediate phase 



between them. Based on the finite size scaling analysis, 
we have shown that the intermediate phase is the criti- 
cal phase with algebraically decaying correlations of the 
order parameter. We found that the phase boundaries 
of the critical phase are of the BKT type. We also ob- 
tained that the numerical values of the critical exponent, 
77, characterizing these two transitions are compatible 
with the theoretical expectations based on the renormal- 
ization group analysis for the six-state clock model<^ It 
should be emphasized that the mapping of the KH model 
onto the six state clock model is valid only below a cer- 
tain temperature which determines the 2D-3D crossover. 
The crossover temperature is, however, larger than the 
temperature of the second BKT transition and, there- 
fore, the low-temperature properties of the KH model 
are captured by the six-state clock model. At high tem- 
peratures, when T is larger than any energy scale in the 
system, the effect of thermal fluctuations is to destroy 
any kind of order and to put the KH model into the 3D 
paramagnetic state. In this regime the mapping of the 
classical KH model to the six state-clock model is not 
valid. 

We also performed a comparative study of honeycomb 
Heisenberg antiferromagnet with the cubic anisotropy 
(|12p . In this case, the continuous symmetry is explicitly 
broken, which allows for the magnetically ordered phase 
to persist up to a finite temperature. We have shown that 
this model has distinct finite temperature properties and 
that the low-T ordered phase is destroyed in a differ- 
ent way. The similarity between the honeycomb Heisen- 
berg antiferromagnet with the cubic anisotropy and the 
KH model is that in both cases the ordered phase is de- 
stroyed in two steps. The main differences between the 
two is that when the cubic anisotropy is taken into ac- 
count explicitly, the intermediate phase is nematic-like, 
and the phase transitions are in the three-states Potts 
and 2D Ising universality classes. 

Finally, we consider the implications of our numer- 
ical results to Na2lr03 and Li2lr03 compounds. We 
have shown that finite temperature magnetic properties 
of these systems can be captured by the extended KH 
model. We found that for the values of the model pa- 
rameters relevant to these systems, the model exhibits 
zigzag magnetic order. Also, our numerical estimates of 
the Neel temperature are close to the experimental val- 
ues. 
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